Introduction

Integrated bulk RNA-seq analysis (MSN9, MSN38, WTC11) of cerebral organoid samples in conditions co-cultured (with/without microglia) and stimulation (ctr/LPS/IFNy).

V1.82: added Line-specific analysis

CO vs COiMG

Organoids co-cultured with vs without microglia.

Variance analysis

PCAplot(pca_cor3_mod1, "COiMg", Shape = "Line", PoV.df=PoV_cor3_mod1, pc.1 = 1, pc.2 = 2) #M30 & M33 are outleirs, let's remove them for now

ggplot(pca.df_upgraded,aes_string(x='COiMG',y='Phagocytosis', fill='COiMG')) +
  stat_compare_means(method = "t.test", label.y = 10, label.x = 1.3) +
  geom_boxplot() +
  geom_boxplot(width=0.1, fill="white") +
  scale_fill_brewer(palette="Dark2") +
  labs(title="Phagocytosis differential PC expression",x="COiMG", y = "PC1 Phagocytosis") +
  theme_classic()

Microglia gene heatmap in ctrl samples

#for control only
ctr.meta <- meta[meta$condition %in% 'ctrl',]
ctr.meta$COiMg[ctr.meta$COiMg %in% 'yes'] <- 'COiMg'
ctr.meta$COiMg[ctr.meta$COiMg %in% 'no'] <- 'CO'

#horizontal
ra <- rowAnnotation(
  COiMG = ctr.meta$COiMg[ctr.meta$condition %in% 'ctrl'],
  col = list(
    COiMG = c("COiMg"='tomato','CO'='lightblue')),
  show_annotation_name = T,
  show_legend = T)


ctr <- t(batch.rem3_mod1[,meta$condition %in% 'ctrl'])
hm_counts <- scale(ctr[,colnames(ctr) %in% gene_panels_subset$Microglia])
ComplexHeatmap::Heatmap(hm_counts,
                        cluster_rows = T,
                        cluster_columns = T,
                        show_row_dend = T,
                        show_column_dend = T,
                        show_row_names = F,
                        show_column_names = F,
                        left_annotation = ra,
                        bottom_annotation = HeatmapAnnotation(
                        text = anno_text(colnames(hm_counts), rot = 300, just = "left",gp=gpar(fontsize=8))),
                        col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(10),
                        name = "Z-score expression")

DEA

Volcano

volcano_plot <- function(res, title = NULL, subtitle = NULL, annotate_by = NULL, type ='ALS', ymax1 = 7.5, ymax2 = 8, xmax1 = -4.5, xmax2 = 4.5){
  res <- 
    mutate(res,
           sig = case_when(
             padj >= 0.05 ~ "non_sig",
             padj < 0.05 & abs(log2FoldChange) < 1 ~ "sig",
             padj < 0.05 & abs(log2FoldChange) >= 1 ~ "sig - strong"
           )) %>%
    mutate(direction = ifelse(log2FoldChange > 0, "up", "down")) %>%
    mutate(log2FoldChange = case_when(
      log2FoldChange > 10 ~ Inf,
      log2FoldChange < -10 ~ -Inf,
      TRUE ~ log2FoldChange
    )) %>%
    mutate(class = paste(sig, direction))
  if( type == "ALS"){
    xpos <- 0.5
    ymax <- ymax1
    xlim <- c(xmax1,xmax2)
  }else{
    xpos <- 0.025
    ymax <- ymax2
    xlim <- c(-0.042, 0.042)
  }
  de_tally <- group_by(res, sig, direction, class) %>% tally() %>%
    filter(sig != "non_sig") %>%
    mutate( position = ifelse(sig == "sig", xpos, 2) ) %>%
    mutate(position = ifelse( direction == "down", -1 * position, position)) %>%
    mutate(n = formatC(n, format="f", big.mark=",", digits=0))
  plot <- res %>%
    mutate( pvalue = ifelse( pvalue < 1e-90, Inf, pvalue)) %>% #threshold at 1e16
    ggplot(aes(x = log2FoldChange, y = -log10(pvalue))) + 
    #geom_point(aes(colour = class ), size = 0.5) +
    rasterise(geom_point(aes(colour = class ), size = 0.5), dpi = 300) +
    scale_colour_manual(values = c("non_sig up" = "gray", 
                                   "non_sig down" = "gray",
                                   "sig up" = "#EB7F56", 
                                   "sig - strong up" = "#B61927",
                                   "sig down" = "#4F8FC4",
                                   "sig - strong down" = "dodgerblue4"
    )) +
    theme_bw() +
    labs(y = expression(-log[10]~P~value), x = expression(log[2]~"(fold change)"), title = title, subtitle = subtitle) +
    guides(colour = FALSE) +
    scale_y_continuous(expand = c(0,0), limits = c(0,ymax)) +
    theme_bw() +
    theme(
      plot.title = element_text(hjust = 0.5),
      plot.subtitle = element_text(hjust = 0.5),
      panel.border = element_blank(),
      axis.ticks = element_line(colour = "black")
    ) +
    geom_text(fontface = "bold", data = de_tally, aes(x = position, y = ymax - 0.5, label = n, colour = class), size = 4) +
    scale_x_continuous(limits = xlim)
  if(!is.null(annotate_by)){
    plot <- plot + 
      ggrepel::geom_text_repel(
        fontface = "italic",
        data = filter(res, symbol %in% annotate_by), 
        aes(x = log2FoldChange, y = -log10(pvalue), label = symbol), 
        min.segment.length = unit(0, "lines"),
        size = 2.3) +
      geom_point(
        data = filter(res, symbol %in% annotate_by), size = 0.8, colour = "black"
      ) +
      geom_point(aes(colour = class ),
                 data = filter(res, symbol %in% annotate_by), size = 0.6
      )
  }
  return(plot)
}
results$COiMG$symbol <- rownames(results$COiMG)
volcano_plot(data.frame(results$COiMG), title = 'CO vs COiMg',
             annotate_by=c('CD36','CD68','CX3CR1','CSF1R','TREM2','TMEM119','ITGAX'), ymax1 = 30, ymax2 = 31)

Top 50 DEGs heatmap

top50.coimg <- rbind(results$COiMG[order(results$COiMG$log2FoldChange, decreasing = T),][results$COiMG[order(results$COiMG$log2FoldChange, decreasing = T),]$padj < 0.05,][1:25,],
results$COiMG[order(results$COiMG$log2FoldChange, decreasing = F),][results$COiMG[order(results$COiMG$log2FoldChange, decreasing = F),]$padj < 0.05,][1:25,])

top50.names <- rownames(top50.coimg)

ra2 <- rowAnnotation(
  COiMG = meta$COiMg,
  col = list(
    COiMG = c("yes"='tomato','no'='lightblue')),
  show_annotation_name = T,
  show_legend = T)


ComplexHeatmap::Heatmap(scale(t(batch.rem3_mod1[rownames(batch.rem3_mod1) %in% top50.names,])),
                        cluster_rows = T,
                        cluster_columns = T,
                        show_row_dend = T,
                        show_column_dend = T,
                        show_row_names = F,
                        show_column_names = F,
                        left_annotation = ra2,
                        bottom_annotation = HeatmapAnnotation(
                        text = anno_text(colnames(t(batch.rem3_mod1[rownames(batch.rem3_mod1) %in% top50.names,])), rot = 300, just = "left",gp=gpar(fontsize=8))),
                        col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(10),
                        name = "Z-score expression")

Which genes of the top DEGs are cilia-associated (CBC)?

setwd('/Users/kubler01/Documents/R projects/gene lists/organoids')
ogenes <- readxl::read_xlsx('gene list_cell type.xlsx', col_names = T)

organoid <- ogenes[,colnames(ogenes) %in% c('CN 1','CN 2','CN 3','CN 4','CN 5','IN','Neuron','Inter','OPC/OL','GPC','NEC 1',
                               'NEC 2','NEC 3','NEC 4','AS1','AS2','AS3','PGC 1','PGC 2','BRC','CBC','UPRC1','UPRC2')][-1,]



coimg.deg <- rbind(results$COiMG[rownames(results$COiMG) %in% na.omit(all_res$COiMG$DEGs[,1]),],
                  results$COiMG[rownames(results$COiMG) %in% na.omit(all_res$COiMG$DEGs[,2]),])


df <- cbind('DEG'=rownames(coimg.deg),'lFC'=round(coimg.deg$log2FoldChange,2),'Cilia-associated'=rownames(coimg.deg) %in% organoid$CBC)
createDT(data.frame(df), "Status", scrollY=1000)

GO

dotplot(all_res$COiMG$GO.up, split="ONTOLOGY", showCategory=10) + facet_grid(ONTOLOGY~., scale="free")

Enrichment

setwd('/Users/kubler01/Documents/R projects/gene lists/organoids')
ogenes <- readxl::read_xlsx('gene list_cell type.xlsx', col_names = T)

organoid <- ogenes[,colnames(ogenes) %in% c('CN 1','CN 2','CN 3','CN 4','CN 5','IN','Neuron','Inter','OPC/OL','GPC','NEC 1',
                               'NEC 2','NEC 3','NEC 4','AS1','AS2','AS3','PGC 1','PGC 2','BRC','CBC','UPRC1','UPRC2')][-1,]


#gene set enrichment analysis:
DEG.enrich.res <- matrix(nrow=2*length(colnames(organoid)), ncol=6)
DEG.rep <- NULL
for (i in names(organoid))
  DEG.rep <- c(DEG.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
DEGs <- list("Up-regulated"=na.omit(all_res$COiMG$DEGs[,1]),"Down-regulated"=na.omit(all_res$COiMG$DEGs[,2]))
rownames(DEG.enrich.res) <- DEG.rep
colnames(DEG.enrich.res ) <- colnames(GSEA.byMod(mod.gl=DEGs, organoid[,i], universe=35324))


for (x in colnames(organoid)){
  DEG.enrich.res[gsub("*_up","",gsub("*_down","",rownames(DEG.enrich.res))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(as.matrix(organoid)[,x]), universe=na.omit(rownames(res))))
}

DEG.enrich.res <- data.frame(DEG.enrich.res)
DEG.enrich.res$log.q <- -log2(DEG.enrich.res[,6])

#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
DEG.enrich <- cbind("Upregulated"=DEG.enrich.res[c(seq(1, length(rownames(DEG.enrich.res))-1, by = 2)),][,7],
                    "Downregulated"=DEG.enrich.res[c(seq(2, length(rownames(DEG.enrich.res)), by = 2)),][,7])
rownames(DEG.enrich) <- names(organoid)

DEG.enrich[DEG.enrich > 20] <- 20
ComplexHeatmap::Heatmap(DEG.enrich,
                        cluster_rows = F,
                        cluster_columns = F,
                        show_row_dend = F,
                        show_column_dend = F,
                        show_row_names = T,
                        col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
                        name = "log2(q)",
)

For microglia pathways.

DEG.enrich.res <- matrix(nrow=2*length(colnames(gene_panels_subset)), ncol=6)
DEG.rep <- NULL
for (i in names(gene_panels_subset))
  DEG.rep <- c(DEG.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
DEGs <- list("Up-regulated"=na.omit(all_res$COiMG$DEGs[,1]),"Down-regulated"=na.omit(all_res$COiMG$DEGs[,2]))
rownames(DEG.enrich.res) <- DEG.rep
colnames(DEG.enrich.res ) <- colnames(GSEA.byMod(mod.gl=DEGs, gene_panels_subset[,i], universe=35324))


for (x in colnames(gene_panels_subset)){
  DEG.enrich.res[gsub("*_up","",gsub("*_down","",rownames(DEG.enrich.res))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(gene_panels_subset[,x]), universe=rownames(res)))
}

DEG.enrich.res <- data.frame(DEG.enrich.res)
DEG.enrich.res$log.q <- -log2(DEG.enrich.res[,6])

#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
DEG.enrich <- cbind("Upregulated"=DEG.enrich.res[c(seq(1, length(rownames(DEG.enrich.res))-1, by = 2)),][,7],
                    "Downregulated"=DEG.enrich.res[c(seq(2, length(rownames(DEG.enrich.res)), by = 2)),][,7])
rownames(DEG.enrich) <- names(gene_panels_subset)
DEG.enrich[DEG.enrich > 20] <- 20
ComplexHeatmap::Heatmap(DEG.enrich,
                        cluster_rows = F,
                        cluster_columns = F,
                        show_row_dend = F,
                        show_column_dend = F,
                        show_row_names = T,
                        col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
                        name = "log2(q)",
)

setwd('/Users/kubler01/Documents/R projects/gene lists/ndd')
cgenes <- readxl::read_xlsx('Gene lists for organoid paper.xlsx', col_names = T)
m <- data.frame(t(toupper(na.omit(cgenes$`ADHD GWAS associated genes through FUMA`))[toupper(na.omit(cgenes$`ADHD GWAS associated genes through FUMA`)) %in% rownames(bulkorg_counts2)]), check.names = F)
for(i in names(cgenes)[-1]){
  k <- na.omit(cgenes[i])
  s <- k[as.matrix(k)[,1] %in% rownames(bulkorg_counts2),]
  m <- rbind.fill(data.frame(m),data.frame(t(s)))
}
rownames(m) <- names(cgenes)
m <- t(m)
#gene set enrichment analysis:
ndd.enrich <- matrix(nrow=2*length(colnames(m)), ncol=6)
ndd.rep <- NULL
for (i in colnames(m))
  ndd.rep <- c(ndd.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
rownames(ndd.enrich) <- ndd.rep
colnames(ndd.enrich) <- colnames(GSEA.byMod(mod.gl=DEGs, m[,1], universe=35324))


for (x in colnames(m)){
  ndd.enrich[gsub("*_up","",gsub("*_down","",rownames(ndd.enrich))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(as.matrix(m)[,x]), universe=na.omit(rownames(res))))
}

ndd.enrich <- data.frame(ndd.enrich)
ndd.enrich$log.q <- -log2(ndd.enrich[,6])

#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
ndd.enrich.ress <- cbind("Upregulated"=ndd.enrich[c(seq(1, length(rownames(ndd.enrich))-1, by = 2)),][,7],
                    "Downregulated"=ndd.enrich[c(seq(2, length(rownames(ndd.enrich)), by = 2)),][,7])
rownames(ndd.enrich.ress) <- colnames(m)

print('No developmental DEGs')
## [1] "No developmental DEGs"

IFNy stimulation

Let’s do this for IFNy stimulation

Variance analysis

PCAplot(pca_cor3_mod1, "condition", Shape = "Line", PoV.df=PoV_cor3_mod1, pc.1 = 1, pc.2 = 2) #M30 & M33 are outleirs, let's remove them for now

ggplot(pca.df_upgraded,aes_string(x='Condition',y='`IFNy`', fill='Condition')) +
  #stat_compare_means(method = "t.test", label.y = 20, label.x = 2) +
  geom_boxplot() +
  geom_boxplot(width=0.1, fill="white") +
  scale_fill_brewer(palette="Dark2") +
  labs(title="Condition differential PC expression",x="Condition", y = "PC1 IFNy response") +
  theme_classic()

DEA

Volcano

results$IFNg$symbol <- rownames(results$IFNg)

volcano_plot(data.frame(results$IFNg), title = 'CTR vs IFNy',
             annotate_by=c('CD36','CD68','CX3CR1','CSF1R','TREM2','TMEM119','ITGAX'), ymax1 = 30, ymax2 = 31)

Top 50 DEGs heatmap

top50.ifng <- rbind(results$IFNg[order(results$IFNg$log2FoldChange, decreasing = T),][results$IFNg[order(results$IFNg$log2FoldChange, decreasing = T),]$padj < 0.05,][1:25,],
results$IFNg[order(results$IFNg$log2FoldChange, decreasing = F),][results$IFNg[order(results$IFNg$log2FoldChange, decreasing = F),]$padj < 0.05,][1:25,])
top50.names2 <- rownames(top50.ifng)
ra3 <- rowAnnotation(
  Stimulation = meta$condition,
  col = list(
    Stimulation = c("ctrl"='lightblue','IFNg'='tomato', 'LPS'='darkred')),
  show_annotation_name = T,
  show_legend = T)


ComplexHeatmap::Heatmap(scale(t(batch.rem3_mod1[rownames(batch.rem3_mod1) %in% top50.names2,])),
                        cluster_rows = T,
                        cluster_columns = T,
                        show_row_dend = T,
                        show_column_dend = T,
                        show_row_names = F,
                        show_column_names = F,
                        left_annotation = ra3,
                        bottom_annotation = HeatmapAnnotation(
                        text = anno_text(colnames(t(batch.rem3_mod1[rownames(batch.rem3_mod1) %in% top50.names2,])), rot = 300, just = "left",gp=gpar(fontsize=8))),
                        col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(10),
                        name = "Z-score expression")

Which genes of the top DEGs are cilia-associated (CBC)?

setwd('/Users/kubler01/Documents/R projects/gene lists/organoids')
ogenes <- readxl::read_xlsx('gene list_cell type.xlsx', col_names = T)

organoid <- ogenes[,colnames(ogenes) %in% c('CN 1','CN 2','CN 3','CN 4','CN 5','IN','Neuron','Inter','OPC/OL','GPC','NEC 1',
                               'NEC 2','NEC 3','NEC 4','AS1','AS2','AS3','PGC 1','PGC 2','BRC','CBC','UPRC1','UPRC2')][-1,]


IFNg.deg <- rbind(results$IFNg[rownames(results$IFNg) %in% na.omit(all_res$IFNg$DEGs[,1]),],
                  results$IFNg[rownames(results$IFNg) %in% na.omit(all_res$IFNg$DEGs[,2]),])


df <- cbind('DEG'=rownames(IFNg.deg),'lFC'=round(IFNg.deg$log2FoldChange,2),'Cilia-associated'=rownames(IFNg.deg) %in% organoid$CBC)
createDT(data.frame(df), "Status", scrollY=1000)

GO

dotplot(all_res$IFNg$GO.up, split="ONTOLOGY", showCategory=10) + facet_grid(ONTOLOGY~., scale="free")

dotplot(all_res$IFNg$GO.down, split="ONTOLOGY", showCategory=10) + facet_grid(ONTOLOGY~., scale="free")

Enrichment

For organoid cell-types.

setwd('/Users/kubler01/Documents/R projects/gene lists/organoids')
ogenes <- readxl::read_xlsx('gene list_cell type.xlsx', col_names = T)

organoid <- ogenes[,colnames(ogenes) %in% c('CN 1','CN 2','CN 3','CN 4','CN 5','IN','Neuron','Inter','OPC/OL','GPC','NEC 1',
                               'NEC 2','NEC 3','NEC 4','AS1','AS2','AS3','PGC 1','PGC 2','BRC','CBC','UPRC1','UPRC2')][-1,]


#gene set enrichment analysis:
DEG.enrich.res <- matrix(nrow=2*length(colnames(organoid)), ncol=6)
DEG.rep <- NULL
for (i in names(organoid))
  DEG.rep <- c(DEG.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
DEGs <- list("Up-regulated"=na.omit(all_res$IFNg$DEGs[,1]),"Down-regulated"=na.omit(all_res$IFNg$DEGs[,2]))
rownames(DEG.enrich.res) <- DEG.rep
colnames(DEG.enrich.res ) <- colnames(GSEA.byMod(mod.gl=DEGs, organoid[,i], universe=35324))


for (x in colnames(organoid)){
  DEG.enrich.res[gsub("*_up","",gsub("*_down","",rownames(DEG.enrich.res))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(as.matrix(organoid)[,x]), universe=na.omit(rownames(res))))
}

DEG.enrich.res <- data.frame(DEG.enrich.res)
DEG.enrich.res$log.q <- -log2(DEG.enrich.res[,6])

#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
DEG.enrich <- cbind("Upregulated"=DEG.enrich.res[c(seq(1, length(rownames(DEG.enrich.res))-1, by = 2)),][,7],
                    "Downregulated"=DEG.enrich.res[c(seq(2, length(rownames(DEG.enrich.res)), by = 2)),][,7])
rownames(DEG.enrich) <- names(organoid)

DEG.enrich[DEG.enrich > 20] <- 20
ComplexHeatmap::Heatmap(DEG.enrich,
                        cluster_rows = F,
                        cluster_columns = F,
                        show_row_dend = F,
                        show_column_dend = F,
                        show_row_names = T,
                        col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
                        name = "log2(q)",
)

For NDD gene lists.

setwd('/Users/kubler01/Documents/R projects/gene lists/ndd')

#gene set enrichment analysis:
ndd.enrich <- matrix(nrow=2*length(colnames(m)), ncol=6)
ndd.rep <- NULL
for (i in colnames(m))
  ndd.rep <- c(ndd.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
rownames(ndd.enrich) <- ndd.rep
colnames(ndd.enrich) <- colnames(GSEA.byMod(mod.gl=DEGs, m[,1], universe=35324))


for (x in colnames(m)){
  ndd.enrich[gsub("*_up","",gsub("*_down","",rownames(ndd.enrich))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(as.matrix(m)[,x]), universe=na.omit(rownames(res))))
}

ndd.enrich <- data.frame(ndd.enrich)
ndd.enrich$log.q <- -log2(ndd.enrich[,6])

#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
ndd.enrich.ress <- cbind("Upregulated"=ndd.enrich[c(seq(1, length(rownames(ndd.enrich))-1, by = 2)),][,7],
                    "Downregulated"=ndd.enrich[c(seq(2, length(rownames(ndd.enrich)), by = 2)),][,7])
rownames(ndd.enrich.ress) <- colnames(m)

ComplexHeatmap::Heatmap(ndd.enrich.ress,
                        cluster_rows = F,
                        cluster_columns = F,
                        show_row_dend = F,
                        show_column_dend = F,
                        show_row_names = T,
                        col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
                        name = "log2(q)",
)

DEGs$`Down-regulated`[DEGs$`Down-regulated` %in% as.character(na.omit(m[,12]))]
## [1] "GRIN2A"

For microglia pathways.

DEG.enrich.res <- matrix(nrow=2*length(colnames(gene_panels_subset)), ncol=6)
DEG.rep <- NULL
for (i in names(gene_panels_subset))
  DEG.rep <- c(DEG.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
DEGs <- list("Up-regulated"=na.omit(all_res$IFNg$DEGs[,1]),"Down-regulated"=na.omit(all_res$IFNg$DEGs[,2]))
rownames(DEG.enrich.res) <- DEG.rep
colnames(DEG.enrich.res ) <- colnames(GSEA.byMod(mod.gl=DEGs, gene_panels_subset[,i], universe=35324))


for (x in colnames(gene_panels_subset)){
  DEG.enrich.res[gsub("*_up","",gsub("*_down","",rownames(DEG.enrich.res))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(gene_panels_subset[,x]), universe=rownames(res)))
}

DEG.enrich.res <- data.frame(DEG.enrich.res)
DEG.enrich.res$log.q <- -log2(DEG.enrich.res[,6])

#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
DEG.enrich <- cbind("Upregulated"=DEG.enrich.res[c(seq(1, length(rownames(DEG.enrich.res))-1, by = 2)),][,7],
                    "Downregulated"=DEG.enrich.res[c(seq(2, length(rownames(DEG.enrich.res)), by = 2)),][,7])
rownames(DEG.enrich) <- names(gene_panels_subset)
DEG.enrich[DEG.enrich > 20] <- 20
ComplexHeatmap::Heatmap(DEG.enrich,
                        cluster_rows = F,
                        cluster_columns = F,
                        show_row_dend = F,
                        show_column_dend = F,
                        show_row_names = T,
                        col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
                        name = "log2(q)",
)

Cell-type marker analysis

Cell-types with only 1 gene found back in our expression data carry principal component expression = 0.

Heatmaps of gene & PC expression

setwd('/Users/kubler01/Documents/R projects/gene lists/organoids')
cell.types <- as.matrix(readxl::read_xlsx('Typical markers cell types organoids 07182023.xlsx', col_names = F))



PC.covs <- matrix(ncol=length(cell.types[,1]),nrow=length(colnames(batch.rem3_mod1)))
colnames(PC.covs) <- cell.types[,1]
for(i in cell.types[,1]){
  pca.cov <- prcomp(t(batch.rem3_mod1[rownames(batch.rem3_mod1) %in% cell.types[cell.types[,1] %in% i,][-1],]))
  PC.covs[,i] <- pca.cov$x[,1]
  #if (table(rownames(batch.rem3_mod1) %in% cell.types[cell.types[,1] %in% i,][-1])[[2]] == 1){
  #  PC.covs[,i] <- t(batch.rem3_mod1[rownames(batch.rem3_mod1) %in% cell.types[cell.types[,1] %in% i,][-1],])
  #}
}
rownames(PC.covs) <- colnames(batch.rem3_mod1)

ra4 <- rowAnnotation(
  Stimulation = meta$condition,
  COiMG = meta$COiMg,
  col = list(COiMG = c('yes'='black','no'='white'),
    Stimulation = c("ctrl"='lightblue','IFNg'='tomato', 'LPS'='darkred')),
  show_annotation_name = T,
  show_legend = T)

ComplexHeatmap::Heatmap(PC.covs,
                        cluster_rows = F,
                        cluster_columns = F,
                        show_row_dend = T,
                        show_column_dend = T,
                        show_row_names = F,
                        show_column_names = F,
                        left_annotation = ra4,
                        bottom_annotation = HeatmapAnnotation(
                        text = anno_text(colnames(PC.covs), rot = 300, just = "left",gp=gpar(fontsize=8))),
                        col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(10),
                        name = "PC1 expression")

#Let's do that for ctr/ifng and coimg/CO only
#for control only
ctr.meta <- meta[meta$condition %in% 'ctrl',]

#horizontal
ctr.ra <- rowAnnotation(
  COiMG = ctr.meta$COiMg,
  col = list(
    COiMG = c("yes"='tomato','no'='lightblue')),
  show_annotation_name = T,
  show_legend = T)


ctr.pca <- PC.covs[meta$condition %in% 'ctrl',]
ComplexHeatmap::Heatmap(ctr.pca,
                        cluster_rows = T,
                        cluster_columns = T,
                        show_row_dend = T,
                        show_column_dend = T,
                        show_row_names = F,
                        show_column_names = F,
                        left_annotation = ctr.ra,
                        bottom_annotation = HeatmapAnnotation(
                        text = anno_text(colnames(ctr.pca), rot = 300, just = "left",gp=gpar(fontsize=8))),
                        col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(10),
                        name = "PC1 expression")

ifng.meta <- meta[!meta$condition %in% 'LPS',]

#horizontal
ifng.ra <- rowAnnotation(
  COiMG = ifng.meta$COiMg,
  Stimulation = ifng.meta$condition,
  col = list(COiMG = c('yes'='black','no'='white'),
    Stimulation = c("ctrl"='lightblue','IFNg'='tomato', 'LPS'='darkred')),
  show_annotation_name = T,
  show_legend = T)


ifng.pca <- PC.covs[!meta$condition %in% 'LPS',]
ComplexHeatmap::Heatmap(ifng.pca,
                        cluster_rows = T,
                        cluster_columns = T,
                        show_row_dend = T,
                        show_column_dend = T,
                        show_row_names = F,
                        show_column_names = F,
                        left_annotation = ifng.ra,
                        bottom_annotation = HeatmapAnnotation(
                        text = anno_text(colnames(ctr.pca), rot = 300, just = "left",gp=gpar(fontsize=8))),
                        col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(10),
                        name = "PC1 expression")

#Let's check only choroid plexus & ependymal
cp.ep <- c(na.omit(cell.types[cell.types[,1] %in% c('Chorioid plexus','Ependymal'),][1,][-1]),
na.omit(cell.types[cell.types[,1] %in% c('Chorioid plexus','Ependymal'),][2,][-1]))

ctr.df <- batch.rem3_mod1[,colnames(batch.rem3_mod1) %in% rownames(ctr.meta)]

ComplexHeatmap::Heatmap(scale(t(ctr.df[rownames(ctr.df) %in% cp.ep,])),
                        cluster_rows = T,
                        cluster_columns = T,
                        show_row_dend = T,
                        show_column_dend = T,
                        show_row_names = F,
                        show_column_names = F,
                        left_annotation = ctr.ra,
                        bottom_annotation = HeatmapAnnotation(
                        text = anno_text(colnames(scale(t(ctr.df[rownames(ctr.df) %in% cp.ep,]))), rot = 300, just = "left",gp=gpar(fontsize=8))),
                        col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(10),
                        name = "Z-score expression")

Boxplots

Let’s make boxplots instead

ifng.pca <- PC.covs[!meta$condition %in% 'LPS',]
ifng.pca_long <- reshape::melt(ifng.pca)
colnames(ifng.pca_long) <- c('Sample','Cell.type','PC')
ifng.pca_long$Stimulation <- ifng.meta$condition
ggplot(ifng.pca_long,aes_string(x='Cell.type',y='PC', fill='Stimulation')) +
  #stat_compare_means(method = "t.test", label.y = 6) +
  geom_boxplot() +
  labs(title="Cell-type marker expression",x="Cell-type", y = "PC1 marker expression") +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, vjust = 0.5))

ctr.meta <- meta[meta$condition %in% 'ctrl',]
ctr.pca <- PC.covs[meta$condition %in% 'ctrl',]
coimg.pca_long <- reshape::melt(ctr.pca)
colnames(coimg.pca_long) <- c('Sample','Cell.type','PC')
coimg.pca_long$Coculture <- ctr.meta$COiMg
ggplot(coimg.pca_long,aes_string(x='Cell.type',y='PC', fill='Coculture')) +
  #stat_compare_means(method = "t.test", label.y = 6) +
  geom_boxplot() +
  labs(title="Cell-type marker expression",x="Cell-type", y = "PC1 marker expression") +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, vjust = 0.5))

Dissecting DEA results

CO-IFNy

These are the results of the DEA of CTR vs IFNy in CO (not co-cultured with microglia) only. Cilia-associated genes remain downregulated. We find this as well in CO vs COiMG and more so also in our single-cell data, which did not have any IFNy-stimulated samples. Meaning: this is a result robust to IFNy AND microglia co-culture separately.

df <- cbind('DEG'=rownames(IFNy.res),'lFC'=round(IFNy.res$log2FoldChange,2),'Cilia-associated'=rownames(IFNy.res) %in% organoid$CBC)
createDT(data.frame(df), "Stimulation", scrollY=1000)

Volcano

new_IFNy$symbol <- rownames(new_IFNy)

volcano_plot(data.frame(new_IFNy), title = 'ctrl vs IFNy',
             annotate_by=c('CD36','CD68','CX3CR1','CSF1R','TREM2','TMEM119','ITGAX'), ymax1 = 30, ymax2 = 31)

CO-COiMg

This is to show that the DEA effect of cilia-associated genes remains present in our bulk-seq analysis when we subset to control samples only and test CO vs COiMg.

df <- cbind('DEG'=rownames(coimg.res),'lFC'=round(coimg.res$log2FoldChange,2),'Cilia-associated'=rownames(coimg.res) %in% organoid$CBC)
createDT(data.frame(df), "COiMG", scrollY=1000)

Volcano

new_coimg$symbol <- rownames(new_coimg)

volcano_plot(data.frame(new_coimg), title = 'CO vs COiMg',
             annotate_by=c('CD36','CD68','CX3CR1','CSF1R','TREM2','TMEM119','ITGAX'), ymax1 = 30, ymax2 = 31)

Correlation of lFCs

overlap <- rownames(coimg.res[rownames(coimg.res) %in% rownames(IFNy.res),])

cor(coimg.res[rownames(coimg.res) %in% overlap,]$log2FoldChange,IFNy.res[rownames(IFNy.res) %in% overlap,]$log2FoldChange)
## [1] 0.718807
createDT(cbind('Genes'=overlap,
               'COiMg'=coimg.res[rownames(coimg.res) %in% overlap,]$log2FoldChange,
               'IFNy'=IFNy.res[rownames(IFNy.res) %in% overlap,]$log2FoldChange), "COiMG", scrollY=1000)

By Line analysis

Correlation of DEA results across lines by condition

IFNy

'IFNy_overlap_wtc11Xmsn9:'
## [1] "IFNy_overlap_wtc11Xmsn9:"
line.res$DEG_IFNy$overlap_wtc11Xmsn9
##  [1] "CPM"      "KIT"      "HPGD"     "JUP"      "KCNJ11"   "AKR1C3"  
##  [7] "LAP3"     "MMP25"    "CD74"     "FAS"      "ARAP2"    "HERPUD1" 
## [13] "PARP12"   "OAS1"     "FKBP5"    "APOL4"    "APOL1"    "TNFSF13B"
## [19] "GIMAP2"   "RIGI"     "ACTA2"    "LGALS3BP" "DHX58"    "CLEC2B"  
## [25] "OAS3"     "OAS2"     "BTN3A3"   "CISH"     "IFIH1"    "GBP3"    
## [31] "IFIT3"    "IFIT2"    "MT2A"     "IFI6"     "APOBEC3F" "BST2"    
## [37] "HELZ2"    "GCH1"     "TRIM22"   "XAF1"     "EPSTI1"   "PLAAT4"  
## [43] "OASL"     "SP110"    "ALDH1L2"  "RTP4"     "IL18BP"   "IFI44L"  
## [49] "IFI44"    "PARP9"    "HERC6"    "PML"      "SECTM1"   "IFITM3"  
## [55] "ITGA9"    "NACC2"    "SLC43A1"  "SLC7A11"  "ANKRD22"  "MOV10"   
## [61] "MX1"      "LY6E"     "GBP2"     "GBP4"     "CDCP1"    "ERAP2"   
## [67] "IFI27"    "BATF2"    "TAP1"     "VAMP5"    "MARCHF3"  "SAMD9L"  
## [73] "PARP10"   "CIITA"    "DDX60L"   "UBA7"     "PYCR1"    "MX2"     
## [79] "USP18"    "SOCS1"    "IRF7"     "IFIT1"    "TRIM69"   "IFITM1"  
## [85] "BTN3A2"   "ISG15"    "RUFY4"    "PLSCR1"   "CASP4"    "HLA-DOA" 
## [91] "HLA-DRA"  "HLA-F"    "IRF9"     "HLA-DPB1" "HLA-DPA1" "HLA-B"   
## [97] "CCDC194"
print('Correlation of IFNy results between WTC11 & MSN9')
## [1] "Correlation of IFNy results between WTC11 & MSN9"
cor(line.res$DEG_IFNy[['WTC11']][rownames(line.res$DEG_IFNy[['WTC11']]) %in% line.res$DEG_IFNy$overlap_wtc11Xmsn9,]$log2FoldChange,
      line.res$DEG_IFNy[['MSN9']][rownames(line.res$DEG_IFNy[['MSN9']]) %in% line.res$DEG_IFNy$overlap_wtc11Xmsn9,]$log2FoldChange)
## [1] 0.2152304
'###################'
## [1] "###################"
'IFNy_overlap_msn9Xmsn38:'
## [1] "IFNy_overlap_msn9Xmsn38:"
line.res$DEG_IFNy$overlap_msn9Xmsn38
##  [1] "MMP25"    "CD74"     "CYP46A1"  "TNS1"     "CMPK2"    "ALDH1L2" 
##  [7] "IL18BP"   "NACC2"    "HPGD"     "MARCHF3"  "CIITA"    "RUFY4"   
## [13] "HLA-DOA"  "HLA-DRA"  "HLA-DPB1" "HLA-DPA1" "LAP3"     "FAS"     
## [19] "ARAP2"    "PARP12"   "DDX3Y"    "SIGLEC1"  "OAS1"     "APOL4"   
## [25] "APOL1"    "TNFSF13B" "GIMAP2"   "RIGI"     "ACTA2"    "LGALS3BP"
## [31] "DHX58"    "CLEC2B"   "OAS3"     "OAS2"     "BTN3A3"   "CISH"    
## [37] "IFIH1"    "GBP3"     "SPP1"     "IFIT3"    "IFIT2"    "MT2A"    
## [43] "ID1"      "IFI6"     "APOBEC3F" "RPS4Y1"   "BST2"     "HELZ2"   
## [49] "GCH1"     "TRIM22"   "XAF1"     "CHI3L1"   "EPSTI1"   "PLAAT4"  
## [55] "OASL"     "SP110"    "RTP4"     "IFI44L"   "IFI44"    "PARP9"   
## [61] "HERC6"    "HERC5"    "PML"      "SECTM1"   "IFITM3"   "ANKRD22" 
## [67] "MOV10"    "MX1"      "LY6E"     "GBP2"     "GBP4"     "CDCP1"   
## [73] "ERAP2"    "IFI27"    "BATF2"    "TAP1"     "VAMP5"    "SAMD9L"  
## [79] "PARP10"   "DDX60L"   "UBA7"     "MX2"      "USP18"    "SOCS1"   
## [85] "IRF7"     "IFIT1"    "TRIM69"   "IFITM1"   "BTN3A2"   "ISG15"   
## [91] "PLSCR1"   "CASP4"    "HLA-F"    "IRF9"     "HLA-B"    "CCDC194"
print('Correlation of IFNy results between MSN38 & MSN9')
## [1] "Correlation of IFNy results between MSN38 & MSN9"
cor(line.res$DEG_IFNy[['MSN38']][rownames(line.res$DEG_IFNy[['MSN38']]) %in% line.res$DEG_IFNy$overlap_msn9Xmsn38,]$log2FoldChange,
    line.res$DEG_IFNy[['MSN9']][rownames(line.res$DEG_IFNy[['MSN9']]) %in% line.res$DEG_IFNy$overlap_msn9Xmsn38,]$log2FoldChange)
## [1] 0.1790641
'###################'
## [1] "###################"
'IFNy_overlap_msn9Xmsn38:'
## [1] "IFNy_overlap_msn9Xmsn38:"
line.res$DEG_IFNy$overlap_msn38Xwtc11
##   [1] "ITM2A"      "HPGD"       "GRIN2A"     "ERVMER34-1" "CFH"       
##   [6] "LAP3"       "CASP10"     "CD38"       "NOS2"       "MMP25"     
##  [11] "IL32"       "ETV7"       "MRC2"       "MVP"        "NUB1"      
##  [16] "WAS"        "CD74"       "SAMD4A"     "TYMP"       "FAS"       
##  [21] "CD44"       "BTN3A1"     "TNFRSF1B"   "TNC"        "LCP2"      
##  [26] "ARAP2"      "GUCA1A"     "EIF2AK2"    "PRDM1"      "PARP12"    
##  [31] "OGFR"       "LZTS1"      "CASP8"      "SBNO2"      "MTHFD2"    
##  [36] "SP100"      "IFI35"      "BCL3"       "ASNS"       "TBX21"     
##  [41] "CACNG5"     "ARHGAP15"   "LAMP3"      "SP140"      "CEACAM1"   
##  [46] "IL12RB2"    "FYB1"       "WDFY1"      "CETP"       "OAS1"      
##  [51] "ICAM1"      "FLT3LG"     "PSME1"      "RNF31"      "JAK2"      
##  [56] "HSD3B7"     "APOL4"      "APOL1"      "PCK2"       "PSME2"     
##  [61] "CD40"       "TRIB3"      "TLDC2"      "SAMHD1"     "SMCHD1"    
##  [66] "TNFSF13B"   "SLC7A5"     "GSDMD"      "ZC3HAV1"    "CAV2"      
##  [71] "CAV1"       "GIMAP2"     "TRIM14"     "RIGI"       "ACTA2"     
##  [76] "MAP3K8"     "LGALS3BP"   "DHX58"      "FAM20A"     "CPZ"       
##  [81] "SLC15A3"    "CARS1"      "CLEC2B"     "OAS3"       "OAS2"      
##  [86] "BTN3A3"     "NCOA7"      "TRIM38"     "TENT5A"     "ITK"       
##  [91] "SLC12A7"    "BCL6"       "CISH"       "OTOF"       "IFIH1"     
##  [96] "STAT1"      "RNF19B"     "FBXO6"      "GBP3"       "GBP1"      
## [101] "IFIT3"      "IFIT2"      "CD274"      "PTK2B"      "TRIM25"    
## [106] "ABCC11"     "TNFSF10"    "NRK"        "NMI"        "BATF3"     
## [111] "DNPEP"      "ZNFX1"      "MT2A"       "IRF1"       "C3"        
## [116] "RBCK1"      "IGFLR1"     "IFI6"       "HELB"       "FGL2"      
## [121] "STEAP4"     "ZFP36"      "APOL3"      "APOL2"      "APOBEC3F"  
## [126] "CHAC1"      "STARD8"     "PLVAP"      "BST2"       "SAG"       
## [131] "HELZ2"      "SHFL"       "IDO1"       "GCH1"       "TRIM21"    
## [136] "TRIM22"     "XAF1"       "EPSTI1"     "PLAAT4"     "DGLUCY"    
## [141] "PSRC1"      "RSAD2"      "IL15RA"     "OASL"       "PRRG4"     
## [146] "PRPH"       "SP110"      "DOCK10"     "ALDH1L2"    "PHF11"     
## [151] "RTP4"       "KLF4"       "RNF144B"    "IL18BP"     "DDX60"     
## [156] "CASP1"      "SQOR"       "IFI44L"     "IFI44"      "STAMBPL1"  
## [161] "PARP9"      "HERC6"      "CXCL9"      "ERP27"      "TAPBPL"    
## [166] "JDP2"       "FBLN5"      "WARS1"      "PML"        "HAPLN3"    
## [171] "MARVELD3"   "NLRC5"      "IRF8"       "SECTM1"     "PMAIP1"    
## [176] "IFITM3"     "MOB3C"      "FCGR2A"     "CYP4V2"     "OSMR"      
## [181] "TACC1"      "NACC2"      "SERPING1"   "AP1S3"      "RASGRP3"   
## [186] "ANKRD22"    "IFIT5"      "LGI4"       "GBP5"       "L3MBTL4"   
## [191] "MOV10"      "PIK3AP1"    "CLIC2"      "ART3"       "MAP3K7CL"  
## [196] "UBE2L6"     "MX1"        "TNFRSF14"   "MEGF11"     "SLAMF8"    
## [201] "IFNAR2"     "CBR3"       "C1R"        "ADAR"       "LY6E"      
## [206] "CXCL16"     "PDPN"       "GBP2"       "GBP4"       "VCAM1"     
## [211] "ATF3"       "CTSS"       "FZD5"       "IFI16"      "CDCP1"     
## [216] "DTX3L"      "LIPH"       "TMEM144"    "IL15"       "ERAP1"     
## [221] "ERAP2"      "TLR3"       "CASP7"      "OTOGL"      "IFI27"     
## [226] "C2"         "PLEKHF1"    "TVP23A"     "B2M"        "TBC1D2B"   
## [231] "NOD2"       "PRRT2"      "BATF2"      "IRF2"       "FILIP1L"   
## [236] "TAP1"       "MLKL"       "STAT3"      "VAMP5"      "LGALS9"    
## [241] "ATF5"       "CXCL10"     "CXCL11"     "PLEKHA2"    "PROKR1"    
## [246] "TRIM56"     "TCAF2"      "STAT2"      "TMEM51"     "ISG20"     
## [251] "CEBPB"      "SLFN11"     "MYD88"      "PARP14"     "PARP15"    
## [256] "C1QB"       "C1QA"       "CD7"        "RNF213"     "MARCHF3"   
## [261] "DES"        "DDIT3"      "A2M"        "SAMD9L"     "PARP10"    
## [266] "EXOC3L1"    "CIITA"      "DDX60L"     "UBA7"       "C1S"       
## [271] "MX2"        "ASCL2"      "FAM162B"    "CSF1"       "SOCS3"     
## [276] "USP18"      "MAFF"       "TNFAIP2"    "SOCS1"      "SP140L"    
## [281] "IRF7"       "IFIT1"      "TRIM69"     "IFITM1"     "BTN3A2"    
## [286] "ISG15"      "EIF4EBP1"   "RUFY4"      "PLSCR1"     "CALHM6"    
## [291] "INSYN2A"    "ANKRD34B"   "HLA-DRB1"   "ARID5A"     "CASP4"     
## [296] "ACSL5"      "GSTK1"      "MAP1LC3C"   "TEAD4"      "TGM2"      
## [301] "HLA-DOA"    "HLA-DMA"    "PSMB8"      "TAP2"       "HLA-DRA"   
## [306] "CARD16"     "LST1"       "HLA-C"      "HLA-E"      "HLA-F"     
## [311] "PSMB10"     "SAMD9"      "ATP10A"     "HLA-A"      "UBD"       
## [316] "IRF9"       "SMTNL1"     "IFI30"      "CEBPD"      "APOL6"     
## [321] "HLA-DPB1"   "C4B"        "PATL2"      "HLA-DPA1"   "CYP21A2"   
## [326] "TAPBP"      "HLA-B"      "KIAA0040"   "OR2I1P"     "APOBEC3G"  
## [331] "PSMB9"      "HLA-DOB"    "HLA-DMB"    "CFB"        "APOBEC3D"  
## [336] "APOBEC3C"   "C4A"        "TMEM158"    "PRSS51"     "TRIL"      
## [341] "ZNF878"     "TRIM34"     "CCDC194"    "CCL5"       "PRAG1"
print('Correlation of IFNy results between MSN38 & WTC11')
## [1] "Correlation of IFNy results between MSN38 & WTC11"
cor(line.res$DEG_IFNy[['MSN38']][rownames(line.res$DEG_IFNy[['MSN38']]) %in% line.res$DEG_IFNy$overlap_msn38Xwtc11,]$log2FoldChange,
    line.res$DEG_IFNy[['WTC11']][rownames(line.res$DEG_IFNy[['WTC11']]) %in% line.res$DEG_IFNy$overlap_msn38Xwtc11,]$log2FoldChange)
## [1] 0.9095827

LinexCondition

for (i in  c('FOXJ1',
             'CLDN5',
             'IDO1',
             'GBP5',
             'MMP10',
             'SOX2',
             'RBFOX3')){
  p <- ggplot(pca.df_upgraded,aes_string(x='LinexCondition',y=i, fill='LinexCondition')) +
    #stat_compare_means(method = "t.test", label.y = 2000, label.x = 1.3) +
    geom_boxplot() +
    geom_boxplot(width=0.1, fill="white") +
    scale_fill_brewer(palette="Dark2") +
    labs(title=paste0("LinexCondition ", i, " expression"),x="LinexCondition", y = paste0(i,' expression')) +
    theme_classic()
  plot(p)
}

COiMg

'IFNy_overlap_wtc11Xmsn9:'
## [1] "IFNy_overlap_wtc11Xmsn9:"
line.res$DEG_COiMg$overlap_wtc11Xmsn9
## [1] "IL1RN"    "IFI44"    "GIPC3"    "ADAMTSL2"
print('Correlation of IFNy results between WTC11 & MSN9')
## [1] "Correlation of IFNy results between WTC11 & MSN9"
cor(line.res$DEG_COiMg[['WTC11']][rownames(line.res$DEG_COiMg[['WTC11']]) %in% line.res$DEG_COiMg$overlap_wtc11Xmsn9,]$log2FoldChange,
      line.res$DEG_COiMg[['MSN9']][rownames(line.res$DEG_COiMg[['MSN9']]) %in% line.res$DEG_COiMg$overlap_wtc11Xmsn9,]$log2FoldChange)
## [1] 0.9463698
'###################'
## [1] "###################"
'IFNy_overlap_msn9Xmsn38:'
## [1] "IFNy_overlap_msn9Xmsn38:"
line.res$DEG_IFNy$overlap_msn9Xmsn38
##  [1] "MMP25"    "CD74"     "CYP46A1"  "TNS1"     "CMPK2"    "ALDH1L2" 
##  [7] "IL18BP"   "NACC2"    "HPGD"     "MARCHF3"  "CIITA"    "RUFY4"   
## [13] "HLA-DOA"  "HLA-DRA"  "HLA-DPB1" "HLA-DPA1" "LAP3"     "FAS"     
## [19] "ARAP2"    "PARP12"   "DDX3Y"    "SIGLEC1"  "OAS1"     "APOL4"   
## [25] "APOL1"    "TNFSF13B" "GIMAP2"   "RIGI"     "ACTA2"    "LGALS3BP"
## [31] "DHX58"    "CLEC2B"   "OAS3"     "OAS2"     "BTN3A3"   "CISH"    
## [37] "IFIH1"    "GBP3"     "SPP1"     "IFIT3"    "IFIT2"    "MT2A"    
## [43] "ID1"      "IFI6"     "APOBEC3F" "RPS4Y1"   "BST2"     "HELZ2"   
## [49] "GCH1"     "TRIM22"   "XAF1"     "CHI3L1"   "EPSTI1"   "PLAAT4"  
## [55] "OASL"     "SP110"    "RTP4"     "IFI44L"   "IFI44"    "PARP9"   
## [61] "HERC6"    "HERC5"    "PML"      "SECTM1"   "IFITM3"   "ANKRD22" 
## [67] "MOV10"    "MX1"      "LY6E"     "GBP2"     "GBP4"     "CDCP1"   
## [73] "ERAP2"    "IFI27"    "BATF2"    "TAP1"     "VAMP5"    "SAMD9L"  
## [79] "PARP10"   "DDX60L"   "UBA7"     "MX2"      "USP18"    "SOCS1"   
## [85] "IRF7"     "IFIT1"    "TRIM69"   "IFITM1"   "BTN3A2"   "ISG15"   
## [91] "PLSCR1"   "CASP4"    "HLA-F"    "IRF9"     "HLA-B"    "CCDC194"
print('Correlation of IFNy results between MSN38 & MSN9')
## [1] "Correlation of IFNy results between MSN38 & MSN9"
cor(line.res$DEG_COiMg[['MSN38']][rownames(line.res$DEG_COiMg[['MSN38']]) %in% line.res$DEG_COiMg$overlap_msn9Xmsn38,]$log2FoldChange,
    line.res$DEG_COiMg[['MSN9']][rownames(line.res$DEG_COiMg[['MSN9']]) %in% line.res$DEG_COiMg$overlap_msn9Xmsn38,]$log2FoldChange)
## [1] 0.6876617
'###################'
## [1] "###################"
'IFNy_overlap_msn9Xmsn38:'
## [1] "IFNy_overlap_msn9Xmsn38:"
line.res$DEG_COiMg$overlap_msn38Xwtc11
##   [1] "BCAS1"    "DBX1"     "FGR"      "TMEM176A" "ITGAL"    "TRAF3IP3"
##   [7] "STAB1"    "CD4"      "BTK"      "MRC2"     "TYROBP"   "ALOX5"   
##  [13] "SLC11A1"  "RUNX3"    "SLAMF7"   "TNFRSF1B" "VCAN"     "MSR1"    
##  [19] "CAPG"     "ADAM28"   "LCP2"     "COL9A2"   "ELN"      "TBXAS1"  
##  [25] "GNA15"    "CD84"     "SPI1"     "APBB1IP"  "EDN1"     "SP140"   
##  [31] "FYB1"     "LAT2"     "NOX4"     "ADAMTS2"  "SIGLEC1"  "RGS1"    
##  [37] "LYZ"      "UNC13D"   "TREM2"    "IL12RB1"  "CYTH4"    "HMOX1"   
##  [43] "NCF4"     "CSF2RB"   "HCK"      "MXRA5"    "ACP5"     "PYCARD"  
##  [49] "AQP9"     "IL4I1"    "RASAL3"   "SIGLEC8"  "CD33"     "PIK3CG"  
##  [55] "TFEC"     "TMEM176B" "ENG"      "ABI3"     "COL1A1"   "IL10RA"  
##  [61] "SLC15A3"  "C11orf21" "BIN2"     "ARHGDIB"  "PTPN6"    "MDFI"    
##  [67] "LOX"      "CD86"     "CYTIP"    "ITGB6"    "PROC"     "PLEK"    
##  [73] "NCF2"     "CD48"     "PADI2"    "SPP1"     "CSF3R"    "TGFBI"   
##  [79] "CCRL2"    "TMIGD3"   "SASH3"    "SRGN"     "BHLHE41"  "ARHGAP9" 
##  [85] "NCKAP1L"  "PMEPA1"   "C3"       "EVI2A"    "ADGRE2"   "FGL2"    
##  [91] "RAC2"     "IRF5"     "WDFY4"    "CD68"     "SIGLEC9"  "APOC1"   
##  [97] "LSP1"     "COL5A1"   "FIBCD1"   "CARD6"    "ALOX5AP"  "CHI3L1"  
## [103] "CHIT1"    "PRAM1"    "GIMAP6"   "GIMAP4"   "ADAMDEC1" "CD180"   
## [109] "PTPN22"   "FST"      "DOCK2"    "SLC37A2"  "HAVCR2"   "SDS"     
## [115] "CD36"     "NT5E"     "DYSF"     "NPL"      "LCP1"     "IL1RN"   
## [121] "TLR4"     "SLCO2B1"  "CASP1"    "ITGA11"   "PLCB2"    "PARVG"   
## [127] "ACVRL1"   "GALNT6"   "GPR65"    "BCL2A1"   "GOLGA6D"  "ITGAX"   
## [133] "PIK3R5"   "VAV1"     "MYO1F"    "CD53"     "FCGR2A"   "LEFTY2"  
## [139] "PTPN7"    "OSBPL10"  "TM4SF19"  "PLA2G7"   "DOK2"     "CDKN2B"  
## [145] "GSN"      "RGS10"    "ST14"     "TAGLN"    "FERMT3"   "FCGR1A"  
## [151] "BMP6"     "SAMSN1"   "SLC7A7"   "VSIG4"    "TDRD9"    "CD109"   
## [157] "SLC34A2"  "TNFRSF14" "NCF1"     "SLAMF8"   "FCER1G"   "C1QC"    
## [163] "RUNX1"    "GOLGA6A"  "ITGB2"    "IL6R"     "LAPTM5"   "FCMR"    
## [169] "HPGDS"    "CTSS"     "S100A11"  "S100A9"   "CCKAR"    "FBLN2"   
## [175] "CCR1"     "TMEM144"  "TAGAP"    "DCSTAMP"  "SYK"      "ALDH1A1" 
## [181] "CYBB"     "HTRA1"    "SPINT1"   "PLD4"     "MFAP4"    "C15orf48"
## [187] "MS4A7"    "MS4A14"   "LAIR1"    "LDLRAD4"  "CXCL8"    "CD52"    
## [193] "ITGAM"    "FPR1"     "GPR34"    "C3AR1"    "SYNPO"    "KLHL6"   
## [199] "C1QB"     "C1QA"     "OLR1"     "TLR10"    "TLR6"     "LRRC25"  
## [205] "NUPR1"    "MSC"      "SSC5D"    "HCLS1"    "CLDN7"    "UBA7"    
## [211] "BGN"      "CSF1R"    "TMEM119"  "APOBR"    "EVI2B"    "CD300LF" 
## [217] "FFAR4"    "ARHGAP30" "TRPV2"    "AMTN"     "DPYD"     "TLR7"    
## [223] "SPN"      "ADAMTSL2" "RCSD1"    "SUCNR1"   "FCGR3A"   "LST1"    
## [229] "STK38L"   "APOC2"    "HLA-DMB"  "CLEC5A"   "CLEC19A"  "PECAM1"  
## [235] "RAB7B"    "CCL3"     "ADORA3"
print('Correlation of IFNy results between MSN38 & WTC11')
## [1] "Correlation of IFNy results between MSN38 & WTC11"
cor(line.res$DEG_COiMg[['MSN38']][rownames(line.res$DEG_COiMg[['MSN38']]) %in% line.res$DEG_COiMg$overlap_msn38Xwtc11,]$log2FoldChange,
    line.res$DEG_COiMg[['WTC11']][rownames(line.res$DEG_COiMg[['WTC11']]) %in% line.res$DEG_COiMg$overlap_msn38Xwtc11,]$log2FoldChange)
## [1] 0.8966225

LinexCOiMG

#LinexCOiMG
for (i in  c('FOXJ1',
             'CLDN5',
             'IDO1',
             'GBP5',
             'MMP10',
             'SOX2',
             'RBFOX3')){
  p <- ggplot(pca.df_upgraded,aes_string(x='LinexCOiMG',y=i, fill='LinexCOiMG')) +
    #stat_compare_means(method = "t.test", label.y = 2000, label.x = 1.3) +
    geom_boxplot() +
    geom_boxplot(width=0.1, fill="white") +
    scale_fill_brewer(palette="Dark2") +
    labs(title=paste0("LinexCOiMG ", i, " expression"),x="LinexCOiMG", y = paste0(i,' expression')) +
    theme_classic()
  plot(p)
}

Interaction analysis

Interaction analysis of co-culture*IFNy

What happens if you co-culture AND stimulate with IFNy?

Variance analysis

pca_cor3_mod1$interaction <- as.factor(pca_cor3_mod1$COiMg):as.factor(pca_cor3_mod1$condition)
PCAplot(pca_cor3_mod1, "interaction", Shape = "Line", PoV.df=PoV_cor3_mod1, pc.1 = 1, pc.2 = 2, 
        colors = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"Spectral")))(6)) #M30 & M33 are outleirs, let's remove them for now

ggplot(pca.df_upgraded,aes_string(x='Condition',y='`IFNy`', fill='Condition')) +
  #stat_compare_means(method = "t.test", label.y = 20, label.x = 2) +
  geom_boxplot() +
  geom_boxplot(width=0.1, fill="white") +
  scale_fill_brewer(palette="Dark2") +
  labs(title="Condition differential PC expression",x="Condition", y = "PC1 IFNy response") +
  theme_classic()

DEA

Volcano

results$COiMGxIFNg$symbol <- rownames(results$COiMGxIFNg)

volcano_plot(data.frame(results$COiMGxIFNg), title = 'CTR vs IFNy',
             annotate_by=c('CD36','CD68','CX3CR1','CSF1R','TREM2','TMEM119','ITGAX'), ymax1 = 10, ymax2 = 11)

DEA boxplots

COiMGxCondition

for (i in  c('FOXJ1',
             'CLDN5',
             'IDO1',
             'GBP5',
             'MMP10',
             'SOX2',
             'RBFOX3')){
  p <- ggplot(pca.df_upgraded,aes_string(x='COiMGxCondition',y=i, fill='COiMGxCondition')) +
    #stat_compare_means(method = "t.test", label.y = 2000, label.x = 1.3) +
    geom_boxplot() +
    geom_boxplot(width=0.1, fill="white") +
    scale_fill_brewer(palette="Dark2") +
    labs(title=paste0("COiMGxCondition ", i, " expression"),x="COiMGxCondition", y = paste0(i,' expression')) +
    theme_classic()
  plot(p)
             }